#include "mex.h"
#include <string.h>
#include<cstdlib>
#include <ctime>
#include <iostream>
using namespace std;
#define INF 1e9
#define MAXSTEPS 100
int findCenter(double* a, double* c, int i, int k)
{
	double min = INF;
	int t = -1;
	for (int j = 0; j < k; j++)
	{
		if ( a[i] - c[j] < min && c[j] - a[i] < min)
		{
			min = a[i] - c[j] > 0 ? a[i] - c[j] : c[j] - a[i];
			t = j;
		}
	}
	return t;
}

void kmeans(double* a, int n, int k, int* b, double* c)
{
	int* p = new int [n];
	memset(p, 0, sizeof(int)*n);	
	// double M = -INF;
	// double m = INF;
	// for(int i = 0; i < n; i++)
	// {
		// if ( a[i] > M ) M = a[i];
		// if ( a[i] < m ) m = a[i];
	// }
	// double h = (M - m) / k;
	// for (int i = 0; i < k; i++)
	// {
		// c[i] = (i + 0.5) * h + m;
	// }
	for (int i = 0; i < k; i++)
	{
		int x = int(n*(double)rand()/RAND_MAX);
		memcpy(c+i,a+x,sizeof(double));
	}
	bool flag = false;
	int step = 0;
	while (!flag)
	{
		int err = 0;
		int* count = new int [k];
		double* s = new double [k];
		memset(count, 0, sizeof(int)*k);
		memset(s, 0, sizeof(double)*k);
		for (int i = 0; i < n; i++)
		{
			int t = findCenter(a, c, i, k); 
			s[t] += a[i];
			count[t]++;
			err += p[i] - t > 0 ? p[i] - t : t - p[i];
			p[i] = t;
		}
		for (int i = 0; i < k; i++)
		{
			if(count[i] == 0)
			{
				for (int i = 0; i < k; i++)
				{
					int x = int(n*(double)rand()/RAND_MAX);
					memcpy(c+i,a+x,sizeof(double));
				}	
				break;
			}
			else
			{
				c[i] = s[i] / count[i];
				if ( err == 0 || step == MAXSTEPS ) flag = true;
			}
		}
		step++;
	}
	memcpy(b, p, sizeof(int)*n);
	delete [] p;
}

double opt(double*a, int n, int k, int*p, double* c)
{
	double s = 0;
	for (int i = 0; i < n; i++)
	{
		s += (a[i]-c[p[i]])*(a[i]-c[p[i]]);
	}
	return s;
}

 //int main()
 //{
	// srand((time(0)));
	// int n = 100;
	// int k = 3;
	// double* a= new double [n];
	// int ITER = 10;
	// for (int i = 0;i < n; i++)
	//	 a[i] = (double)rand()/RAND_MAX;
	// double min = INF;
	// int* tmpb = new int [n];
	// double* tmpc = new double [k];
	// int* b = new int [n];
	// double* c = new double [k];
	//
	// for (int i = 0 ; i < ITER; i++)
	// {	
	//	 kmeans(a, n, k, tmpb, tmpc);
	//	 double optVal = opt(a, n, k, tmpb, tmpc);
	//	 if ( optVal < min)
	//	 {
	//		 min = optVal;
	//		 memcpy(b, tmpb, sizeof(int)*n);
	//		 memcpy(c, tmpc, sizeof(double)*k);
	//	 }
	// }
	// system("pause");
	// return 0;
 //}


void mexFunction(int nlhs, mxArray* plhs[], int nrhs, mxArray* prhs[])
{
	if ( nrhs < 2)
		mexErrMsgTxt("2 input arguments is required.");
	else if (nrhs > 3)
		mexErrMsgTxt("Too many input arguments.");
	else if( nlhs > 2)
		mexErrMsgTxt("Too many output arguments.");
	
	int n = mxGetM(prhs[0])*mxGetN(prhs[0]);
	int k = (int)mxGetScalar(prhs[1]);
	plhs[0] = mxCreateNumericMatrix(1,k, mxDOUBLE_CLASS, mxREAL);
	plhs[1] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[0]), 
		mxGetDimensions(prhs[0]), mxINT32_CLASS, mxREAL);
	double* a = (double*)mxGetPr(prhs[0]);
	int ITER = 10;
	if ( nrhs == 3)
		ITER = mxGetScalar(prhs[2]);
	if ( ITER < 1 )
		mexErrMsgTxt("Iteration number must be larger than 1.");

	double min = INF;
	int* tmpb = new int [n];
	double* tmpc = new double [k];
	int* b = (int *)mxGetPr(plhs[1]);
	double* c = mxGetPr(plhs[0]);
	for (int i = 0 ; i < ITER; i++)
	{	
		kmeans(a, n, k, tmpb, tmpc);
		double optVal = opt(a, n, k, tmpb, tmpc);
		if ( optVal < min)
		{
			min = optVal;
			memcpy(b, tmpb, sizeof(int)*n);
			memcpy(c, tmpc, sizeof(double)*k);
		}
	}
}






























